if (!requireNamespace("BiocManager", quietly = TRUE))
  install.packages("BiocManager")
library(BiocManager)

library(GEOquery)
library(limma)
library(ggplot2)
library(aplot)
library(ggrepel)

setwd("D:\\ProgramFiles\\R\\Work\\mRNA")

source("CIBERSORT.R")

load("var.RData")

set.seed(2022)

#miRNA是算好数的文件直接下载后本地读取。用getGEO获取一下分组信息

fileName <- dir() 
mydata <- read.table(fileName[1], header = T, sep = "\t", row.names = 1)
length(unique(rownames(mydata)))
ges <- data.frame(row.names = rownames(mydata))


# 读取其余文件，并逐一并到第一个文件上
for(i in 1:(length(fileName)-1)){
   data <- read.csv2(fileName[i],header = T, sep = "\t", row.names = 1)
   ges[,substr(fileName[i],start = 1, stop = 10)] <- data[,"Counts"]
}


# 读取自己下载的信息，由于没有对照，所以没能成功
# var_exp <- read.delim("GSE92696_non-normalized.txt.gz", stringsAsFactors = F)
# var_exp <- read.ilmn("GSE92696_non-normalized.txt.gz", expr = "SAMPLE", 
#                      probeid = "ID_REF", other.columns = "Detection Pval")

# var <- neqc(var_exp, detection.p = "Detection Pval")
# lumi包对lumina芯片数据的读取和标准化质控
# library(lumi)
# x <- lumiR.batch("GSE92696_non-normalized.txt.gz")
# data.eset <- lumiExpresso(x.lumi) #质控
# data.exprs <- exprs(data.eset)

ges <- getGEO("GSE22165", GSEMatrix = T)
ges <- getGEO(filename = "GSE92696_series_matrix.txt.gz", GSEMatrix = T, getGPL = F)

exper <- read.delim("GSE92696_series_matrix.txt.gz", comment.char = "!", header = T)
annot <- read.delim("GPL10558.annot.gz", skip = "28", header = T)
annot2 <- data.frame(annot$Gene.symbol, annot$ID)
group <- read.delim("GSE92696_series_matrix.txt.gz", skip = 28, nrow = 10, header = T)
group <- group[,-1]
groups <- as.factor(ifelse(group[9,] %in% "neurological outcome (cerebral performance category): CPC1",0,1))


colnames(annot2) <- c("symbol", "ID_REF")
mx <- merge(exper, annot2, all.x = T, sort = T,  )


mx <- as.data.frame(exprs(ges[[1]]))
mx$symbol <- featureData(ges[[1]])[["Gene Symbol"]]

#根据Symbol去重，选择在所有样本中表达量之和最大的探针。
library(dplyr)
mx1 <- mx %>% mutate(rowMean = rowMeans(.[grep("GSM", names(.))])) %>%
  arrange(desc(rowMean)) %>%
  distinct(symbol, .keep_all = TRUE) %>%
  select(-rowMean)

rownames(mx1) <- mx1$symbol
mx2 <- select(mx1, -symbol)
mx2 <- select(mx2, -ID_REF)

my <- mx2[fina_gene,]
my <- mx2[inge,]

groups <-phenoData(ges[[1]])$description.1
groups <- factor(rep(c("GOOD", "BAD"), each = 5), levels = c ("GOOD", "BAD"))
#去掉小于3个样本表达的探针
ges[ges==0] <- NA
mx2 <- ges[rowSums(!is.na(ges)) > 3,]

sum(is.na(mx2))

#KNN法补全缺失值
library(DMwR2)
my <- knnImputation(mx2)

#上个代码报错，采用多重插补法补全缺失值
library(mice)
impdatas <- mice(mx2, seed = 2022)
my <- complete(impdatas,1)
my <- data.matrix(my)

sum(is.na(my))

tmy <- data.frame(t(my))
groups <- as.factor(ifelse(groups %in% c("CPC: 1", "CPC: 2"),"GOOD","BAD"))
groups <- factor(rep(c("GOOD", "BAD"), each = 25), levels = c ("GOOD", "BAD"))

tmy$groups <- groups

#PCA分析
pca<- prcomp(tmy,scale. = T, retx = T)
plot(pca$x[,1],pca$x[,2])

#注意stat_ellipse绘制椭圆，type = "norm"假定元素呈多元正态分布
df_pca <- data.frame(pca$x, groups = groups)
pdf("pca.pdf",height = 6, width = 8)
ggplot(df_pca,aes(x = PC1, y = PC2, color = groups)) +
  geom_point(size = 3)+
  theme_classic()+
  theme(legend.position = "none")+
  stat_ellipse(type = "norm", 
               alpha = 0.3, geom = "polygon", aes(fill = groups),
               linetype = "blank", size = 0.5)+
  #geom_label(data = df_pca, aes(label = rownames(df_pca)))+
  xlab(paste0("PC1 (",round(pca$sdev[1]/sum(pca$sdev)*100,2),"%)"))+
  ylab(paste0("PC2 (",round(pca$sdev[2]/sum(pca$sdev)*100,2),"%)"))
  #annotate('text', label = 'BAD', x = 155, y = 120, size = 6, colour = '#00bfc4') +
  #annotate('text', label = 'GOOD', x = -200, y = 120, size = 6, colour = '#f8766d') 
dev.off()

#去掉偏离过远的异常值。
myp<-myp[,!colnames(myp) %in% c("GSM731335", "GSM731371")]
which(colnames(mx2) == "GSM731371")
groups <- groups[-61]

#去除异常值后重新做差异分析，limma法，适用于芯片数据
my <- my[,!colnames(my) %in% c("GSM731335", "GSM731371")]

#背景矫正和标准化，有的需要有的不需要。
my <- backgroundCorrect.matrix(mx2,method = "normexp")
my <- normalizeBetweenArrays(my, method = "quantile")

#有的需要取log，有的不需要。
my2 <- log2(my)

#绘制箱型图看表达情况
par(cex = 0.7)
n.sample = ncol(my)
if (n.sample>40) par(cex = 0.5)
cols <-rainbow(n.sample*1.2)
boxplot(my, col=cols, main = "expression value", las=2)

#设置分组矩阵
design <- model.matrix(~0+groups)
colnames(design) = levels(groups)
rownames(design) = colnames(my)

#构建比较矩阵
contrast.matrix <-makeContrasts(paste0(unique(groups),collapse = "-"),
                                 levels = design)

#差异分析
fit <- lmFit(my,design)
fit2 <- contrasts.fit(fit,contrast.matrix)
fit2 <- eBayes(fit2,0.01)
tempOutput = topTable(fit2,adjust = "fdr", coef = 1, n = Inf)
nrDEG = na.omit(tempOutput)
head(nrDEG)
sum(nrDEG$adj.P.Val<0.05)
sum(nrDEG$logFC>0.5)
#绘制火山图
nrDEG$Threshold<-as.factor(ifelse(nrDEG$adj.P.Val<0.05&
                                    abs(nrDEG$logFC)>=0.5,
                                        ifelse(nrDEG$logFC>0.5,"UP","DOWN"),"NO"))
range(nrDEG$logFC)
#nrDEG2<-nrDEG[-1,]
nrDEG$ID <- rownames(nrDEG)
pdf("hsp.pdf", height = 10, width = 12)
ggplot(data = nrDEG, aes(x=logFC, 
                         y=-log10(adj.P.Val),colour = Threshold ))+
  geom_point(alpha=0.5,size = 4)+
  scale_color_manual(values = c("blue","grey","red"))+
  theme_bw() + xlim(-2,2)+
  theme_classic()+
  theme(
    axis.title = element_text(size = 15),
    axis.text = element_text(size = 12),
    legend.title = element_text(size = 15),
    legend.text = element_text(size = 12),
    legend.position = c(1,1),
    legend.justification = c(1,1),
    legend.background = element_rect(fill = NULL, colour = "black"))+
  geom_vline(xintercept = c(-0.5,0.5),lty = 2, col = "black", lwd = 0.4)+
  geom_hline(yintercept = -log10(0.05),lty = 2, col = "black", lwd = 0.4)+
  geom_text_repel(data = subset(nrDEG, adj.P.Val<0.05&abs(nrDEG$logFC) > 0.5),
                  aes(label = ID))
dev.off()

#使用DESeq2对高通量数据进行差异分析，使用未补全的数据mx2或者使用补全缺失值的my，
#最终选择了还是使用补全了缺失值的方法。
library(DESeq2)
mx2[is.na(mx2)] <- 0
df_grous <- data.frame(colnames(mx2),groups)
dds <- DESeqDataSetFromMatrix(countData = my, colData = df_grous, design = ~groups )
dds <- DESeq(dds)
res <- results(dds, contrast = c("groups", "GOOD","BAD"))
res$padj[is.na(res$padj)] <-1 
sum(res$padj < 0.05)
sum(res$log2FoldChange > 1)

#绘制火山图
df_hsd <- data.frame(res$log2FoldChange, res$padj, row.names = rownames(res))
df_hsd$ID <- rownames(df_hsd)
df_hsd$Threshold<-factor(ifelse(df_hsd$res.padj < 0.05&
                                    abs(df_hsd$res.log2FoldChange) > 0.5,
                                   ifelse(df_hsd$res.log2FoldChange> 0.5,"UP","DOWN"),"NO"))
pdf("hsp.pdf", height = 10, width = 12)
ggplot(df_hsd,aes(x= res.log2FoldChange, y = -log10(res.padj), color = Threshold))+
  geom_point(alpha=0.5,size = 4)+
  scale_color_manual(values = c("blue","grey","red"))+
  theme_bw() + xlim(-4,4)+
  theme_classic()+
  theme(
    axis.title = element_text(size = 15),
    axis.text = element_text(size = 12),
    legend.title = element_text(size = 15),
    legend.text = element_text(size = 12),
    legend.position = c(1,1),
    legend.justification = c(1,1),
    legend.background = element_rect(fill = NULL, colour = "black"))+
  geom_vline(xintercept = c(-0.5,0.5),lty = 2, col = "black", lwd = 0.4)+
  geom_hline(yintercept = -log10(0.05),lty = 2, col = "black", lwd = 0.4)+
  geom_text_repel(data = subset(df_hsd, res.padj < 0.05&
                                  abs(res.log2FoldChange) > 0.5),
                  aes(label = ID))+
  labs(x = "logFC", y = "-log10(adj.P.Val)")
dev.off()

#制作绘制热图所需数据表
pID<-nrDEG[nrDEG$adj.P.Val<0.05,]$ID
uID<-cbind(pID)
rownames(uID)<-pID
heatmx<-merge(my2,uID,by="row.names")
hdf <-as.data.frame(heatmx)
hdf <- hdf[-36,]
library(reshape2)
hdfc <-melt(hdf,id=c("pID","Row.names"))
hdfc[is.na(hdfc)]<-0 

#绘制热图,unclass把row.names数据类型解开
hp <- ggplot(hdfc, aes(variable, unclass(hdfc[,2])))+
  geom_tile(aes(fill=value))+
  scale_fill_gradient(low="white",high ='red')+
  #热图变点图
  #geom_point(aes(size=value,color=value))+
  #scale_color_gradient(low = "white", high = "red")+
  theme_classic(base_family = 'serif')+
  theme(axis.text.x = element_text(angle = -45,hjust = 0),
        axis.title.x = element_blank(),
        panel.border = element_rect(fill = NA, color = "black", size = 0.5, linetype = "solid"),
        axis.line = element_line(colour = "black",size = 1)
        )+
  geom_vline(xintercept = 5.5, size = 0.5)+
  ylab('Circulating microRNAs')


#绘制分组条形图
dfg<-data.frame(groups = groups, sampale = colnames(my2),y = rep(c("gr"),times =10))
dfgp<-ggplot(dfg,aes(x= sampale, y = y,fill = groups))+
  geom_tile() +
  theme_void()+
  labs(fill = "groups")+
  scale_fill_discrete(labels = c("CPC 1-2","CPC 3-5"),breaks = c("good","bad"))+
  theme(legend.title = element_text(family = "serif"),
        legend.text = element_text(family = "serif"))

#使用Y叔的aplot将图拼起来
library(aplot)
pdf("heatmap_miRNA.pdf",height = 8, width = 8)
insert_top(hp,dfgp, height=0.05)
dev.off()


#碎石图
#载荷因子princomp用loadings，只适用于行大于列的数据，prcomp用rotation。
screeplot(pca, type = "lines")
pca_pr <- predict(pca)
summary(pca_pr, loadings = TRUE)
loadings(pca_pr)

#选择累积方差大于80%的作为主成分
i = 1
while (sum(pca$sdev[1:i]/sum(pca$sdev))<0.8) {
  i = i+1
  print(i)
  }
#取前90个PC
sum(pca$sdev[1:90]/sum(pca$sdev))
pca_u <- pca$rotation[,1:90]
#——————————————————分割线————————————————#
#没有用提取的大于80%的PC。因为无法选取差异基因，只是利用PCA删除了异常的样本，
#再次进行差异分析，并用差异表达的基因，重新绘制得到了PCA图。
#下一步利用差异基因进行WGCNA分析。
#BiocManager::install("WGCNA")
#library(WGCNA)
#用MEGENA进行基因共表达分析

library(MEGENA)
#删除变量，只剩一个my
rm(list=setdiff(ls(),c("my", "groups")))
# data(Sample_Expression)
# input parameters
n.cores <- 16; # number of cores/threads to call for PCP
doPar <-TRUE; # do we want to parallelize?
method = "pearson" # method for correlation. either pearson or spearman. 
FDR.cutoff = 0.05 # FDR threshold to define significant correlations upon shuffling samples. 
module.pval = 0.05 # module significance p-value. Recommended is 0.05. 
hub.pval = 0.05 # connectivity significance p-value based random tetrahedral networks
cor.perm = 10; # number of permutations for calculating FDRs for all correlation pairs. 
hub.perm = 100; # number of permutations for calculating connectivity significance p-value. 

# annotation to be done on the downstream
annot.table=NULL
id.col = 1
symbol.col= 2


#### register multiple cores if needed: note that set.parallel.backend() is deprecated. 
run.par = doPar & (getDoParWorkers() == 1)
if (run.par)
{
  cl <- parallel::makeCluster(n.cores)
  registerDoParallel(cl)
  # check how many workers are there
  cat(paste("number of cores to use:",getDoParWorkers(),"\n",sep = ""))
}

ijw <- calculate.correlation(my,doPerm = cor.perm,output.corTable = FALSE,output.permFDR = FALSE)
start.time <- Sys.time()
##### calculate PFN
el <- calculate.PFN(ijw[,1:3],doPar = doPar,num.cores = n.cores,keep.track = FALSE)
end.time <- Sys.time()
time.taken <- end.time - start.time
time.taken

g <- graph.data.frame(el,directed = FALSE)

##### perform MCA clustering.
MEGENA.output <- do.MEGENA(g,
                           mod.pval = 0.05, hub.pval = 0.05,remove.unsig = TRUE,
                           min.size = 10,max.size = vcount(g)/2,
                           doPar = doPar,num.cores = n.cores,n.perm = hub.perm,
                           save.output = FALSE)

###### unregister cores as these are not needed anymore.
if (getDoParWorkers() > 1)
{
  env <- foreach:::.foreachGlobals
  rm(list=ls(name=env), pos=env)
}

summary.output <- MEGENA.ModuleSummary(MEGENA.output,
                                       mod.pvalue = 0.05,hub.pvalue = 0.05,
                                       min.size = 10,max.size = vcount(g)/2,
                                       annot.table = NULL,id.col = 1,symbol.col = 2,
                                       output.sig = TRUE)

if (!is.null(annot.table))
{
  # update annotation to map to gene symbols
  V(g)$name <- paste(annot.table[[symbol.col]][match(V(g)$name,annot.table[[id.col]])],V(g)$name,sep = "|")
  summary.output <- output[c("mapped.modules","module.table")]
  names(summary.output)[1] <- "modules"
}

print(summary.output$module.table)

pnet.obj <- plot_module(output.summary = summary.output,PFN = g,subset.module = "c1_2",
                        layout = "kamada.kawai",label.hubs.only = TRUE,
                        gene.set = NULL,color.code =  "grey",
                        output.plot = FALSE,out.dir = "modulePlot",col.names = c("magenta","green","cyan"),label.scaleFactor = 20,
                        hubLabel.col = "black",hubLabel.sizeProp = 1,show.topn.hubs = Inf,show.legend = TRUE)

print(pnet.obj[[1]])

module.table <- summary.output$module.table
colnames(module.table)[1] <- "id" # first column of module table must be labelled as "id".

hierarchy.obj <- plot_module_hierarchy(module.table = module.table,label.scaleFactor = 0.15,
                                       arrow.size = 0.03,node.label.color = "blue")

print(hierarchy.obj[[1]])

#hub节点信息，我们所要的关键靶点就是这个了。
hub.summary <- get.hub.summary(MEGENA.output)

#所有输入基因的计算结果，这个里面的is.hub应该是小范围内的，不是全局那种。
module.df <- module_convert_to_table(MEGENA.output, mod.pval = 0.05, hub.pval = 0.05,
                                     min.size = 10, max.size = 5000)

#差异分析的数据和MEGENA的数据取交集
DEGs <- load("DEGs.RData")
rm(list = setdiff(ls(), c ("nrDEG", "hub.summary","tmy","groups")))
inge <- intersect(subset(nrDEG, nrDEG$P.Value<0.1)$ID, hub.summary$node)
inge <- intersect(subset(nrDEG, nrDEG$P.Value<0.05)$ID, module.df$id)
inge <- intersect(subset(df_hsd, df_hsd$res.padj<0.05)$ID, hub.summary$node)

#miRNA中DEGs与MEGENA的交集只有三个，所以舍弃MEGENA的结果进行下面分析。
inge <- subset(df_hsd, df_hsd$res.padj < 0.05)$ID
inge <- gsub('[-]', '.', inge)
input <- tmy[, (colnames(tmy) %in% inge)]
input <- cbind(groups, input)

dim(input)
#机器学习,svm-RFE。并行代码。行为样本，列为基因，第一列为分组。
library(e1071)
library(Rmpi)
library(snow)
library(parallel)

source("D:\\ProgramFiles\\R\\Work\\msvmRFE.R")

nfold = 10
nrows = nrow(input)
folds = rep(1:nfold, len=nrows)[sample(nrows)]
folds = lapply(1:nfold, function(x) which(folds == x))

#make a cluster
cl <- makeMPIcluster(mpi.universe.size())

clusterExport(cl, list("input","svmRFE","getWeights","svm"))
results <-parLapply(cl,folds, svmRFE.wrap, input, k=10, halve.above=100)
top.features = WriteFeatures(results, input, save=F)

clusterExport(cl, list("top.features","results", "tune","tune.control"))
featsweep = parLapply(cl,1:100, FeatSweep.wrap, results, input)
stopCluster(cl)

no.info = min(prop.table(table(input[,1])))
errors = sapply(featsweep, function(x) ifelse(is.null(x), NA, x$error))
#dev.new(width=4, height=4, bg='red')
pdf("svm_rfe.pdf", height = 8, width = 10)
PlotErrors(errors, no.info=no.info)
dev.off()
plot(top.features)
mpi.exit()


#LASSO回归，输入只能是矩阵。行为样本，列为基因，分组只能是（0,1），多组可以转化MultinomialExample
library(glmnet)
# data("BinomialExample")
# x <- BinomialExample$x
# y <- BinomialExample$y
x <- input[,-1]
levels(groups) <- c(0, 1)
y <- as.integer(as.character(groups))
fit <- glmnet(x, y, family = "binomial", nlambda = 100, alpha = 1)
plot(fit, xvar = "lambda", label = TRUE)
cvfit <- cv.glmnet(data.matrix(x), y, nfolds = 10)
cvfit$lambda.min
cvfit$lambda.1se

pdf("lasso.pdf", height = 8, width = 10)
plot(cvfit)
dev.off()
coef <- coef(cvfit, s = "lambda.min")
factors <- as.matrix(coef)
factors[factors == 0] <- NA
factors <- na.omit(factors)
write.csv(factors, "lasso.csv")

#随机森林，行为样本，列为基因，搞一列为分类 
library(randomForest)

infr <- as.data.frame(input)
colnames(infr) <- make.names(colnames(infr))
infr[,1] <- as.factor(infr[,1])
model <- randomForest(groups~., data = infr, proximity = TRUE, ntree = 1000)

#画图看种多少树合适
oob.error.data <- data.frame(
  Trees=rep(1:nrow(model$err.rate), times=3),
  Type=rep(c("OOB", "GOOD", "BAD"), each=nrow(model$err.rate)),
  Error=c(model$err.rate[,"OOB"],
          model$err.rate[,"GOOD"],
          model$err.rate[,"BAD"]))
ggplot(data=oob.error.data, aes(x=Trees, y=Error)) +
  geom_line(aes(color=Type))

#确定mtry值
oob.values <- vector(length=20)
for(i in 1:20) {
  temp.model <- randomForest(groups ~ ., data = infr, mtry=i, ntree=1000)
  oob.values[i] <- temp.model$err.rate[nrow(temp.model$err.rate),1]
}
oob.values

#最终确定的结果
model <- randomForest(as.factor(groups)~.,data = infr, mtry = 3, ntree = 600,
                      proximity = TRUE, importance = TRUE)
rfi1 <- importance(model, type = 1)
rfi2 <- importance(model, type = 2)

rfi1 <- rfi1[order(rfi1, decreasing = T),]
rfi2 <- rfi2[order(rfi2, decreasing = T),]

rfif <- intersect(names(rfi1[1:30]), names(rfi2[1:30]))
pdf ("randf.pdf", height = 8, width = 10)
varImpPlot(model)
dev.off()

fina_gene <- intersect(top.features[1:68,]$FeatureName, rownames(factors))
fina_gene <- intersect(fina_gene, rfif)
fina_gene

#绘制韦恩图
library(venn)
venn_list <- list(top.features[1:68,]$FeatureName, rownames(factors), rfif)
pdf("venn.pdf", height = 4, width = 5)
venn(venn_list, snames = c("SVM-REF", "LASSO", "RF"), zcolor = c( "#F8766D", "#00BA38", "#619CCF"),
     ellipse = T, box = F, col = c( "#F8766D", "#00BA38", "#619CCF"), ilcs = 1, sncs = 1)
dev.off()

#基因名中的.和-互换。
fina_gene <- gsub('[.]', '-', fina_gene)
#因为其他基因名无-5p所以可以粗暴的去掉。
fina_gene <- fina_gene<- c("hsa-let-7a" , "hsa-let-7c",  "hsa-miR-3138",   "hsa-miR-545-5p" ,"hsa-miR-483-5p")
#看验证数据集中是否包含fina_gene
sum(grepl("313",rownames(mx2)))
mx2[grep("let",rownames(mx2)),]

tmy <- cbind(tmy, mx2["hsa-let-7a",])

#计算var数据集中final基因的表达在两组间的比较
wilcox <- wilcox.test(tmy[,"PRKX"]~group,  data = tmy)
gewilcox <- function(x){
  test <- wilcox.test(tmy[,x]~group, data = tmy)
  pval <- test$p.value
  return(pval)
}
gewilcox("CFLAR.PRKX")

tmy$CFLAR.PRKX <- tmy$CFLAR - tmy$PRKX

sapply(fina_gene, gewilcox)
#画柱状图
expr_df <- data.frame(PRKX = tmy$PRKX, CFLAR = tmy$CFLAR, groups = tmy$group, 
                      CFLAR.PRKX = tmy$CFLAR.PRKX)
expr_df <- melt(expr_df, id = "groups", variable.name = "gene", value.name = "expr")
expr_df$gene <- factor(expr_df$gene, levels = c("CFLAR", "PRKX", "CFLAR.PRKX"))

expr_seg <- aggregate(expr_df$expr, by = list(expr_df$groups, expr_df$gene), max) %>%
  pivot_wider(names_from = Group.1, values_from = x)
seg_sign <- c(1,3)
expr_seg<- expr_seg[seg_sign,]
expr_seg$xx <- seg_sign
colnames(expr_seg) <- c("gene", "GOOD", "BAD", "xx")

#上下各画一半然后拼起来即可。
library(aplot)
pdown <-ggplot(expr_df, aes(x = gene, y = expr))+
  geom_boxplot(aes(fill = groups))+
  theme_classic()+
  theme(legend.title=element_blank())+
  scale_fill_discrete(labels = c("GOOD","BAD"))+
  scale_x_discrete(labels = c("CFLAR", "PRKX", "CFLAR-PRKX"))+
  labs(x = "", y = "Expression")+
  geom_segment(data = expr_seg, aes(x = xx-0.2,xend =xx +0.2, y = ifelse(GOOD > BAD, GOOD, BAD)+0.3,
                                   yend = ifelse(GOOD > BAD, GOOD, BAD)+0.3),
               lineend = "round")+
  geom_segment(data = expr_seg, aes(x = xx+0.2,xend =xx +0.2, y = BAD +0.2,
                                   yend = ifelse(GOOD > BAD, GOOD, BAD)+0.3))+
  geom_segment(data = expr_seg, aes(x = xx-0.2,xend =xx -0.2, y = GOOD+0.2,
                                   yend = ifelse(GOOD > BAD, GOOD, BAD)+0.3))+
  annotate("text", x = expr_seg$xx ,y = ifelse(expr_seg$GOOD > expr_seg$BAD,
                                              expr_seg$GOOD, expr_seg$BAD)+0.3,
           label = "*",size = 8, color = "red")+
  ylim(c(5,7.5))
pup <- ggplot(expr_df, aes(x = gene, y = expr))+
  geom_boxplot(aes(fill = groups))+
  theme_classic()+
  theme(legend.title=element_blank(), axis.text.x = element_blank(), 
        axis.ticks.x = element_blank(), axis.line.x = element_blank())+
  scale_fill_discrete(labels = c("GOOD","BAD"))+
  scale_x_discrete(labels = c("CFLAR", "PRKX", "CFLAR-PRKX"))+
  labs(x = "", y = "")+
  geom_segment(data = expr_seg, aes(x = xx-0.2,xend =xx +0.2, y = ifelse(GOOD > BAD, GOOD, BAD)+0.3,
                                    yend = ifelse(GOOD > BAD, GOOD, BAD)+0.3),
               lineend = "round")+
  geom_segment(data = expr_seg, aes(x = xx+0.2,xend =xx +0.2, y = BAD +0.2,
                                    yend = ifelse(GOOD > BAD, GOOD, BAD)+0.3))+
  geom_segment(data = expr_seg, aes(x = xx-0.2,xend =xx -0.2, y = GOOD+0.2,
                                    yend = ifelse(GOOD > BAD, GOOD, BAD)+0.3))+
  annotate("text", x = expr_seg$xx ,y = ifelse(expr_seg$GOOD > expr_seg$BAD,
                                               expr_seg$GOOD, expr_seg$BAD)+0.3,
           label = "*",size = 8, color = "red")+
  ylim(c(12,13.5))

pdf("expr.pdf", height = 6, width = 6)
insert_top(pdown,pup,height = 0.6)
dev.off()

#绘制ROC曲线
rm(list = setdiff(ls(), "fina_gene"))
library(pROC)
rocobj <- plot.roc(tmy$groups, tmy[,2], auc = T)
text(0.25, 0.4, labels = paste("miR-483-5p_AUC = ", round(rocobj$auc, 3)),
     col = "black")

aggregate(tmy[,1:7], list(tmy$groups),mean)

#计算每一个基因的AUC面积。
TP <- rep(0,10)
for (i in 1:2){
  rocobj <- roc(tmy$groups, tmy[,i], smooth = T)
  TP[i] <- rocobj$auc
}
max(TP)
which.max(TP)

#AUC最好的是基因中表达增加AUC最大的减去表达减小的AUC最大的。并绘图。
pdf("ROC.pdf", height = 6, width = 6)
rocobj <- plot.roc(tmy$groups, tmy[,5], main = "ROC",
                   smooth = F, col = "#F8766D", xlim = c(1,0),
                   levels=c("GOOD", "BAD"),direction="<" )
rocobj1 <- lines.roc(tmy$groups, tmy[,2], smooth = F,
                     col = "#619CCF", levels=c("GOOD", "BAD"),direction="<" )
rocobj2 <- lines.roc(tmy$groups, tmy[,1], smooth = F,
                     col = "#00BA38", levels=c("GOOD", "BAD"),direction="<" )
text(0.25, 0.4, labels = paste("miR-483_AUC = ", round(rocobj$auc, 3)),
     col = "#F8766D")
text(0.25, 0.35, labels = paste("let-7a_AUC = ", round(rocobj2$auc, 3)),
     col = "#00BA38")
text(0.25, 0.3, labels = paste("let-7c_AUC = ", round(rocobj1$auc, 3)),
     col = "#619CCF")
legend("bottomright", legend=c("miR-483", "let-7a", "let-7c"),
       col=c( "#F8766D", "#00BA38", "#619CCF"), lwd=2)
dev.off()

#免疫浸润分析CIBERSORT，数据格式为第一列为基因名+Gene symbol，不能有空值和重复，
#重复基因取平均值或最大值处理，第一行为样本名称，数据为data.frame。不能去log

write.table(my, "my.txt", sep = "\t", row.names = T, col.names = T)
myjrr <- CIBERSORT("LM22.txt", "my.txt", perm = 1000, QN = T)

#构建画图数据框
library(tidyr)
df_my <- data.frame(myjrr[,1:22]) %>% mutate(groups = groups) %>%
  mutate(sample = rownames(myjrr)) %>% pivot_longer(cols = colnames(.)[1:22],
                                                    names_to = "cell_type",
                                                    values_to = 'value')
#构建排序，按表达量排
plot_order <- df_my[df_my$groups=="GOOD",] %>% 
  group_by(cell_type) %>% 
  summarise(m = median(value)) %>% 
  arrange(desc(m)) %>% 
  pull(cell_type)

df_my$cell_type <- factor(df_my$cell_type,levels = plot_order)
df_my$groups <- factor(df_my$groups, levels = c("GOOD", "BAD"))
#构造线段数据框
df_seg <- aggregate(df_my$value, by = list(df_my$groups, df_my$cell_type), max) %>%
  pivot_wider(names_from = Group.1, values_from = x)
seg_sign <- c(1,3,7,11,12,18,19,21)
df_segg <- df_seg[seg_sign,]
df_segg$xx <- seg_sign


#画图
library(ggpubr)
pdf("myjr.pdf", height = 9, width = 12)
ggplot(df_my, aes(cell_type, value))+
  geom_boxplot(aes(fill = groups))+
  theme_classic()+
  theme(legend.position = "top",legend.title=element_blank(),
        axis.text.x = element_text(angle = 45, hjust = 1 ),)+
  labs(y = "Cell composition", x = NULL)+
  # stat_compare_means(aes(group =  groups),#这个是ggpubr的加星标。下面自己添加的。
  #                    label = "p.signif",
  #                    method = "wilcox.test",
  #                    hide.ns = T,size = 5)+
  geom_segment(data = df_segg, aes(x = xx-0.2,xend =xx +0.2, y = ifelse(GOOD > BAD, GOOD, BAD)+0.03,
                                   yend = ifelse(GOOD > BAD, GOOD, BAD)+0.03),
               lineend = "round")+
  geom_segment(data = df_segg, aes(x = xx+0.2,xend =xx +0.2, y = BAD +0.02,
                                   yend = ifelse(GOOD > BAD, GOOD, BAD)+0.03))+
  geom_segment(data = df_segg, aes(x = xx-0.2,xend =xx -0.2, y = GOOD+0.02,
                                   yend = ifelse(GOOD > BAD, GOOD, BAD)+0.03))+
  annotate("text", x = df_segg$xx ,y = ifelse(df_segg$GOOD > df_segg$BAD,
                                              df_segg$GOOD, df_segg$BAD)+0.03,
           label = "*",size = 8, color = "red")
dev.off()  

#单基因与免疫浸润细胞的关系
cor_fun <- function(y){
  x = "PRKX"
  cor_x <- my[x,]
  cor_y <- myjrr[,y]
  if(sd(cor_y) == 0){cor_y[1] = 0.0000001}
  cor_t <- cor.test(cor_x, cor_y, method = "pearson")
  r = cor_t[["estimate"]][["cor"]]
  p = cor_t[["p.value"]]
  return(c(x,y,r,p))
}
cor_fun("Monocytes")
prkx_my <- sapply(colnames(myjrr), cor_fun)

#绘制棒棒糖图
df_bbt <- data.frame(t(cflar_my))
df_bbt$X3 <- as.numeric(df_bbt$X3)
df_bbt$X4 <- as.numeric(df_bbt$X4)
df_bbt$X2 <- reorder(df_bbt$X2, -df_bbt$X3)
pdf("cflar_bb.pdf", height = 6, width = 9)
ggplot(df_bbt, aes(x = X3, y = X2))+
  geom_segment(aes(x = 0, xend = X3, y = X2, yend = X2), col = ifelse(df_bbt$X4 < 0.05, "red", "black"))+
  geom_point(aes(size = abs(X3),col = ifelse(df_bbt$X4 < 0.05, "p < 0.05", "p > 0.05")))+
  labs(x = "Correlation Coefficient", y = NULL, size = "abs(r)", color = "pvalue")+
  theme_classic()+
  geom_vline(xintercept = c(-0.3,0.3),linetype = "dashed", lwd = 0.4, colour = "grey50")
dev.off()
#使用person相关计算与关键基因表达有关的基因
library(compiler) #compfun()可以把自己写的函数编译一下，运行更快
library(future.apply)#future.apply可以多线程
#配合plan(multisession)
cor_df <- c()
cor_fun <- function(y){
   x = "CFLAR"
   data = my
   cor_x <- data[x,]
   cor_y <- data[y,]
   cor_t <- cor.test(cor_x, cor_y, method = "pearson")
   r = cor_t[["estimate"]][["cor"]]
   p = cor_t[["p.value"]]
   ifelse(abs(r) > 0.6 & p <0.05,cor_df <- cbind(x, y, r, p), cor_df <- cbind(NA, NA, NA, NA))
   return(cor_df)
}

cor_df <- cor_fun("NPL")

cor_df <- sapply(rownames(my), cor_fun)

cor_dff <- na.omit(t(cor_df))
cor_dfc <- na.omit(t(cor_df))
write.table(cor_dfc, "CFLAR.txt", col.names = F, row.names = F, quote = F, sep = "\t")

#根据共表达和miRNA的结果做GSEA分析，（没结果）
library(clusterProfiler)
library(org.Hs.eg.db)
library(ReactomePA)
library(msigdbr)
data(geneList, package = "DOSE")
anage <- read.csv("Merged Network_2 default node.csv", header = T)
listge <- merge(anage, nrDEG, by.x = "name", by.y = "ID")
geneid <- bitr(listge$name, fromType = "SYMBOL", toType = "ENTREZID", OrgDb = "org.Hs.eg.db")
geneid <- rbind(geneid,c("QARS", 5859))
genefc <- merge(listge,geneid, by.x = "name", by.y = "SYMBOL")

list <- genefc$logFC
names(list) <- genefc$ENTREZID
list <- sort(list, decreasing = T)
#开始富集
go <- gseGO(list, OrgDb = "org.Hs.eg.db", minGSSize = 5,ont = "ALL", pvalueCutoff = 0.2, verbose = FALSE)
kk <- gseKEGG(list, OrgDb = "org.Hs.eg.db" )
re <- gsePathway(list,minGSSize = 1, pvalueCutoff = 1, pAdjustMethod = "BH"
                 , verbose = FALSE, readable = T)
dotplot(go)
go <- enrichGO(genefc$ENTREZID, ont = "ALL", OrgDb = org.Hs.eg.db)


msigdbr_species()
Hs_msigdbr <- msigdbr(species="Homo sapiens")
colnames(Hs_msigdbr)
Hs_df <- as.data.frame(Hs_msigdbr[,c('gs_name','entrez_gene','gene_symbol')])
head(Hs_df)
em_msig <- GSEA(list,TERM2GENE=Hs_df[,c(1,2)], pvalueCutoff = 1,by="DOSE", nPerm = 10)
dotplot(em_msig)
wk <- gseWP(list, organism = "Homo sapiens", by="DOSE")

write.csv(re,file = "re.csv")
